1 AR Models

${\displaystyle y_t}$ is regressed on its own lagged values ${\displaystyle y_{t-1},\cdots ,y_{t-p}}$.

For predicting ${\displaystyle y_{n+1}}$, plug ${\displaystyle t=n+1}$ into (1.1): $$y_{n+1}={\hat{\varphi }}_0+{\hat{\varphi }}_1{y}_n+\cdots +{\hat{\varphi }}_p{y}_{n+1-p}.$$
Note that ${\displaystyle y_n,y_{n-1},\cdots ,y_{n+1-p}}$ are all observed and they are the last ${\displaystyle p}$ observations. Then for ${\displaystyle y_{n+2}}$, $$y_{n+2}={\hat{\varphi }}_0+{\hat{\varphi }}_1{y}_{n+1}+\cdots +{\hat{\varphi }}_p{y}_{n+2-p}.$$
Here ${\displaystyle y_{n+1}}$ is not observed, but we can replace with predicted value ${\displaystyle {\hat{y}}_{n+1}}$. And recursively, predict value for any time point.

We can see AR as simply regression of the observed time series on lagged versions of itself.

2 Relation to MLE

Start with ${\displaystyle \mathrm{AR}(1)}$: $$\begin{array}{}\text{(2.1)}& y_t={\varphi }_0+{\varphi }_1{y}_{t-1}+{\epsilon }_t,t=2,\cdots ,n.\end{array}$$

2.1 Usual Regression

(2) looks just like $$\begin{array}{}\text{(2.2)}& y_i={\beta }_0+{\beta }_1{x}_i+{\epsilon }_i.\end{array}$$ For this, given independence, $$\begin{array}{rl}\mathrm{Likelihood}& =f_{x_1,y_1,\cdots ,x_n,y_n|\theta }(x_1,y_1,\cdots ,x_n,y_n)\\ & =\prod _{i=1}^n{f}_{x_i,y_i|\theta }(x_i,y_i)=\prod _{i=1}^n{f}_{y_i|x_i,\theta }(y_i)f_{x_i|\theta }(x_i).\end{array}$$
Now replace ${\displaystyle y_i}$ by ${\displaystyle {\beta }_0+{\beta }_1{x}_i+{\epsilon }_i}$: $$\begin{array}{rl}\mathrm{Likelihood}& =\prod _{i=1}^n{f}_{{\beta }_0+{\beta }_1{x}_i+{\epsilon }_i|x_i,\theta }(y_i)f_{x_i|\theta }(x_i)\\ & =\prod _{i=1}^n{f}_{{\epsilon }_i|x_i,\theta }(y_i-{\beta }_0-{\beta }_1{x}_i)f_{x_i|\theta }(x_i)\\ & =\prod _{i=1}^m{f}_{{\epsilon }_i|\theta }(y_i-{\beta }_0-{\beta }_1{x}_i)f_{x_i|\theta }(x_i)\\ & =\prod _{i=1}^m\frac{1}{\sqrt{2\pi }\sigma }\exp (-\frac{(y_i-{\beta }_0-{\beta }_1{x}_i)}{2{\sigma }^2})f_{x_i|\theta }(x_i)\\ & ={\left(\frac{1}{\sqrt{2\pi }\sigma }\right)}^m\exp (-\frac{1}{2{\sigma }^2}\sum _{i=1}^m(y_i-{\beta }_0-{\beta }_1{x}_i{)}^2)\prod _{i=1}^m{f}_{x_i|\theta }(x_i)\\ & \propto {\left(\frac{1}{\sqrt{2\pi }\sigma }\right)}^m\exp (-\frac{1}{2{\sigma }^2}\sum _{i=1}^m(y_i-{\beta }_0-{\beta }_1{x}_i{)}^2).\end{array}$$
To sum up, here are the assumptions we made:

• Independence of ${\displaystyle (x_1,y_1),\cdots ,(x_n,y_n)}$.
• Model equation (linearity)
• Independence of ${\displaystyle {\epsilon }_i,x_i}$
• ${\displaystyle {\epsilon }_i\sim \mathcal{N}(0,{\sigma }^2)}$
• Density of ${\displaystyle x_i}$ independent of ${\displaystyle \theta =({\beta }_0,{\beta }_1,\sigma )}$.

2.2 Back to ${\displaystyle {\displaystyle \mathrm{AR}(1)}}$

When we try to apply (3) to (2), these assumptions don't hold:

• Independence of ${\displaystyle (x_i,y_i)}$
• Density of ${\displaystyle x_t=y_{t-1}}$ depends on ${\displaystyle \theta }$ (${\displaystyle {\varphi }_0,{\varphi }_1,\sigma }$ affects ${\displaystyle y_{t-1}}$)

So now given ${\displaystyle y_1,\cdots ,y_n}$, $$\begin{array}{rl}\text{Likelihood for (2)}& =f_{y_1,\cdots ,y_n|\theta }(y_1,\cdots ,y_n)\\ & =f_{y_1|\theta }(y_1)f_{y_2|y_1,\theta }(y_2)\cdots f_{y_n|y_1,\cdots ,y_{n-1},\theta }(y_n)\\ & =f_{y_1|\theta }(y_1)\prod _{t=2}^n{f}_{y_t|y_1,\cdots ,y_{t-1},\theta }(y_t)\\ & =f_{y_1|\theta }(y_1)\prod _{t=2}^n{f}_{{\epsilon }_t|y_1,\cdots ,y_{t-1},\theta }(y_t-{\varphi }_0-{\varphi }_1{y}_{t-1})\\ & =f_{y_1|\theta }(y_1)\prod _{t=2}^n{f}_{{\epsilon }_t}(y_t-{\varphi }_0-{\varphi }_1{y}_{t-1})\\ & =f_{y_1|\theta }(y_1)\prod _{t=2}^n\frac{1}{\sqrt{2\pi }\sigma }\exp (-\frac{1}{2{\sigma }^2}(y_t-{\varphi }_0-{\varphi }_1{y}_{t-1}{)}^2)\\ & =f_{y_1|\theta }(y_1){\left(\frac{1}{\sqrt{2\pi }\sigma }\right)}^{n-1}\exp (-\frac{1}{2{\sigma }^2}\sum _{t=2}^n(y_t-{\varphi }_0-{\varphi }_1{y}_{t-1}{)}^2).\end{array}$$
Here we assume

• Model equation (2.1).
• Independence of ${\displaystyle {\epsilon }_t}$ and ${\displaystyle y_1,\cdots ,y_{t-1}}$ for ${\displaystyle t=2,\cdots ,n-1}$.
• ${\displaystyle {\epsilon }_t\sim \mathcal{N}(0,{\sigma }^2)}$.

2.2.1 Computation of ${\displaystyle {\displaystyle f_{y_1|\theta }(y_1)}}$

We can't get ${\displaystyle f_{y_1|\theta }(y_1)}$ from equation. There are two approaches to get it.

First, we simply assume ${\displaystyle f_{y_1|\theta }(y_1)}$ does not depend on ${\displaystyle \theta }$. So $$\begin{array}{}\text{(2.3)}& \mathrm{Likelihood}\propto {\left(\frac{1}{\sqrt{2\pi }\sigma }\right)}^{n-1}\exp (-\frac{1}{2{\sigma }^2}\sum _{t=2}^n(y_t-{\varphi }_0-{\varphi }_1{y}_{t-1}{)}^2).\end{array}$$
It's easy to verify that the estimates ${\displaystyle {\hat{\varphi }}_0,{\hat{\varphi }}_1}$ are identical to those obtained in (2.1). We call this Conditional MLE.

Second, we extend the model to ${\displaystyle t=1,0,-1,-2,\cdots }$. So $$\begin{array}{rl}{y}_1& ={\varphi }_0+{\varphi }_1{y}_0+{\epsilon }_1\\ & ={\varphi }_0+{\varphi }_1({\varphi }_0+{\varphi }_1{y}_{-1}+{\epsilon }_0)+{\epsilon }_1\\ & ={\varphi }_0(1+{\varphi }_1)+{\varphi }_1^2{y}_{-1}+{\varphi }_1{\epsilon }_0+{\epsilon }_1\\ & ={\varphi }_0(1+{\varphi }_1)+{\varphi }_1^2({\varphi }_0+{\varphi }_1{y}_{-2}+{\epsilon }_{-1})+{\varphi }_1{\epsilon }_0+{\epsilon }_1\\ & ={\varphi }_0(1+{\varphi }_1+{\varphi }_1^2)+{\varphi }_1^3{y}_{-2}+{\varphi }_1^2{\epsilon }_{-1}+{\varphi }_1{\epsilon }_0+{\epsilon }_1\\ & =\cdots ={\varphi }_0\sum _{j=0}^M{\varphi }_1^j+{\varphi }_1^{M+1}{y}_{-M}+\sum _{j=0}^M{\varphi }_1^j{\epsilon }_{1-j}.\end{array}$$
If ${\displaystyle |{\varphi }_1|<1}$, coefficient ${\displaystyle {\varphi }_1^{M+1}}$ is very small, so $$y_1\approx {\varphi }_0\sum _{j=0}^M{\varphi }_1^j+\sum _{j=0}^M{\varphi }_1^j{\epsilon }_{1-j}\approx {\varphi }_0\sum _{j=0}^{\mathrm{\infty }}{\varphi }_1^j+\sum _{j=0}^{\mathrm{\infty }}{\varphi }_1^j{\epsilon }_{1-j}=\frac{{\varphi }_0}{1-{\varphi }_1}+\sum _{j=0}^{\mathrm{\infty }}{\varphi }_1^j{\epsilon }_{1-j}.$$
Since ${\displaystyle \mathbb{E}\left(\sum _{j=0}^{\mathrm{\infty }}{\varphi }_1^j{\epsilon }_{1-j}\right)=0}$, and $$\mathrm{Var}\left(\sum _{j=0}^{\mathrm{\infty }}{\varphi }_1^j{\epsilon }_{1-j}\right)=\sum _{j=0}^{\mathrm{\infty }}\mathrm{Var}({\varphi }_1^j{\epsilon }_{1-j})=\sum _{j=0}^{\mathrm{\infty }}{\varphi }_1^{2j}\mathrm{Var}({\epsilon }_{1-j})={\sigma }^2\sum _{j=0}^{\mathrm{\infty }}{\varphi }_1^{2j}=\frac{{\sigma }^2}{1-{\varphi }_1^2}.$$
Thus when ${\displaystyle |{\varphi }_1|<1}$, $$y_1\sim \mathcal{N}(\frac{{\varphi }_0}{1-{\varphi }_1},\frac{{\sigma }^2}{1-{\varphi }_1^2}),$$ so $$f_{y_1|\theta }(y_1)=\frac{\sqrt{1-{\varphi }_1^2}}{\sqrt{2\pi }\sigma }\exp (-\frac{1-{\varphi }_1^2}{2{\sigma }^2}{(y_1-\frac{{\varphi }_0}{1-{\varphi }_1})}^2),$$ so finally $$\begin{array}{}\text{(5)}& \mathrm{Likelihood}=\frac{\sqrt{1-{\varphi }_1^2}}{(\sqrt{2\pi }\sigma {)}^n}\exp (-\frac{1-{\varphi }_1^2}{2{\sigma }^2}{(y_1-\frac{{\varphi }_0}{1-{\varphi }_1})}^2-\frac{1}{2{\sigma }^2}\sum _{t=2}^n(y_t-{\varphi }_0-{\varphi }_1{y}_{t-1}{)}^2).\end{array}$$
Now compare (4) and (5). (5) is referred to as full likelihood for ${\displaystyle \mathrm{AR}(1)}$ and therefore full MLE. (4) and (5) will be quite close when ${\displaystyle |{\varphi }_1|<1}$ and ${\displaystyle n}$ is large.

2.3 ${\displaystyle {\displaystyle \mathrm{AR}(p)}}$

The ${\displaystyle \mathrm{AR}(p)}$ model is given by $$y_t={\varphi }_0+{\varphi }_1{y}_{t-1}+\cdots +{\varphi }_p{y}_{t-p}+{\epsilon }_t.$$
The likelihood is $$f_{y_1,\cdots ,y_n|\theta }(y_1,\cdots ,y_n)=f_{y_{p+1},\cdots ,y_n|y_1,\cdots ,y_p,\theta }(y_{p+1},\cdots ,y_n)f_{y_1,\cdots ,y_p|\theta }(y_1,\cdots ,y_p).$$
The conditional likelihood is $$\begin{array}{rl}& f_{y_{p+1},\cdots ,y_n|y_1,\cdots ,y_p,\theta }(y_{p+1},\cdots ,y_n)=\prod _{t=p+1}^n{f}_{y_t|y_{t-1},\cdots ,y_1}(y_t)\\ =& \prod _{t=p+1}^n{f}_{{\epsilon }_t|y_{t-1},\cdots ,y_1}(y_t-{\varphi }_0-{\varphi }_1{y}_{t-1}-\cdots -{\varphi }_p{y}_{t-p})\\ =& {\left(\frac{1}{\sqrt{2\pi }\sigma }\right)}^{n-p}\exp (-\frac{1}{2{\sigma }^2}\sum _{t=p+1}^n(y_t-{\varphi }_0-{\varphi }_1{y}_{t-1}-\cdots -{\varphi }_p{y}_{t-p}{)}^2).\end{array}$$
Here we assume ${\displaystyle {\epsilon }_t|y_{t-1},\cdots ,y_1\sim \mathcal{N}(0,{\sigma }^2),t=p+1,\cdots ,n}$.

To obtain the parameter estimates, we can directly maximize the likelihood. And since ${\displaystyle f_{y_1,\cdots ,y_p|\theta }(y_1,\cdots ,y_p)}$ does not depend on ${\displaystyle \theta }$, it is equivalent to maximizing the conditional likelihood.

2.3.1 Bayesian Approach

However, if we want to derive ${\displaystyle f_{y_1,\cdots ,y_p|\theta }(y_1,\cdots ,y_p)}$ in a more principled way, we have to use (1.1) for smaller values of ${\displaystyle t}$ but it's complicated and not really worth it. We can also work under some "stationarity" assumptions on ${\displaystyle {\varphi }_0,\cdots ,{\varphi }_p}$ (much simpler than conditional likelihood): use matrix notation (see here), and $$\mathrm{likelihood}{\propto }_{\theta }{\left(\frac{1}{\sqrt{2\pi }\sigma }\right)}^{n-1}\exp (-\frac{||Y-X\beta |{|}^2}{2{\sigma }^2}).$$ Assume ${\displaystyle {\varphi }_0,\cdots ,{\varphi }_p,\log \sigma \stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}}{\sim }\mathrm{Unif}(-C,C)}$, then by here, $$\beta |\mathrm{data}\sim t_{n-2p-1,p+1}(\hat{\beta },{\hat{\sigma }}^2(X^{\mathrm{T}}X{)}^{-1}),$$ where ${\displaystyle \hat{\beta }=(X^{\mathrm{T}}X{)}^{-1}{X}^{\mathrm{T}}Y,\hat{\sigma }=\sqrt{\frac{||Y-X\hat{\beta }|{|}^2}{n-2p-1}}}$. If inference for ${\displaystyle \sigma }$ is desired, we can use $$\frac{||Y-X\hat{\beta }|{|}^2}{{\sigma }^2}|\mathrm{data}\sim {\chi }_{n-2p-1}^2.$$

Bayesian inference for AR is identical to linear regression models because of the same likelihood. Bayesian inference only cares about the likelihood.
Frequentist inference is based on MLE, given by ${\displaystyle \hat{\beta }}$ and ${\displaystyle {\hat{\sigma }}_{\mathrm{MLE}}=\sqrt{\frac{||Y-X\hat{\beta }|{|}^2}{n-p}}}$. The analysis is quite different from linear regression. The results are slightly different but close.

3 Predictions and Difference Equations

Given a fitted ${\displaystyle \mathrm{AR}(p)}$ model with ${\displaystyle {\hat{\varphi }}_0,\cdots ,{\hat{\varphi }}_p}$, predictions ${\displaystyle {\hat{y}}_{n+i}}$ for ${\displaystyle i=1,2,\cdots }$ are obtained by: $$\begin{array}{}\text{(3.1)}& {\hat{y}}_{n+i}={\hat{\varphi }}_0+{\hat{\varphi }}_1{\hat{y}}_{n+i-1}+\cdots +{\hat{\varphi }}_p{\hat{y}}_{n+i-p},i=1,2,\cdots \end{array}$$ where the recursion is initialized with $${\hat{y}}_j=y_j,j=n,n-1,\cdots ,n+1-p.$$
Or we can rewrite to $$\begin{array}{}\text{(3.2)}& u_k={\alpha }_0+{\alpha }_1{u}_{k-1}+\cdots +{\alpha }_p{u}_{k-p},k=p,p+1,\cdots \end{array}$$ (Initialized by ${\displaystyle u_0,\cdots ,u_{p-1}}$). (3.2) is called a difference equation of order ${\displaystyle p}$.

3.1 First Order (${\displaystyle {\displaystyle p=1}}$)

Now (3.2) becomes $$u_k={\alpha }_0+{\alpha }_1{u}_{k-1}$$ along with initialized ${\displaystyle u_0}$. Convert to a homogeneous equation (with no intercept term) by taking ${\displaystyle v_k=u_k-\frac{{\alpha }_0}{1-{\alpha }_1}}$: $$v_k={\alpha }_1{v}_{k-1}.$$
So $$v_k={\alpha }_1^k{v}_0\Rightarrow u_k=\frac{1-{\alpha }_1^k}{1-{\alpha }_1}{\alpha }_0+{\alpha }_1^k{u}_0.$$

• ${\displaystyle |{\alpha }_1|<1}$: ${\displaystyle u_k}$ converges exponentially to ${\displaystyle \frac{{\alpha }_0}{1-{\alpha }_1}}$;
• ${\displaystyle |{\alpha }_1|>1}$: explodes to infinity exponentially;
• ${\displaystyle {\alpha }_1=1}$: ${\displaystyle u_k=k{\alpha }_0+u_0}$;
• ${\displaystyle {\alpha }_1=-1}$: ${\displaystyle u_k}$ oscillates between ${\displaystyle u_0}$ and ${\displaystyle {\alpha }_0-u_0}$.

3.2 General Case: Bayesian Approach

In the Bayesian context, prediction is done via joint probability distribution of ${\displaystyle y_{n+1},\cdots ,y_{n+k}}$ conditional on ${\displaystyle y_1,\cdots ,y_n}$. Now consider conditional expectations: $$\begin{array}{}\text{(3.3)}& \mathbb{E}(y_{n+i}|y_1,\cdots ,y_n)=\int \mathbb{E}(y_{n+i}|y_1,\cdots ,y_n,\theta )f_{\theta |y_1,\cdots ,y_n}(\theta )\mathrm{d}\theta .\end{array}$$
First calculate ${\displaystyle {\hat{y}}_{n+i}(\theta )=\mathbb{E}(y_{n+i}|y_1,\cdots ,y_n,\theta )}$ for fixed ${\displaystyle \theta }$: $$\begin{array}{}\text{(3.4)}& {\hat{y}}_{n+i}(\theta )={\varphi }_0+{\varphi }_1{\hat{y}}_{n+i-1}(\theta )+\cdots +{\varphi }_p{\hat{y}}_{n+i-p}(\theta ).\end{array}$$
If we initialize this with ${\displaystyle {\hat{y}}_j(\theta )=y_j}$, ${\displaystyle j=n,n-1,\cdots ,n+1-p}$, then (3.4) can be evaluated in sequence for ${\displaystyle i=1,2,\cdots }$.

Now (3.3) becomes $$\mathbb{E}(y_{n+i}|y_1,\cdots ,y_n)=\int {\hat{y}}_{n+i}(\theta )f_{\theta |y_1,\cdots ,y_n}(\theta )\mathrm{d}\theta .$$
We can do one of two things to compute:

1. Generate posterior samples ${\displaystyle {\theta }^{(1)},\cdots ,{\theta }^{(N)}}$ from ${\displaystyle f_{\theta |y_1,\cdots ,y_n}(\theta )}$, then $$\mathbb{E}(y_{n+i}|y_1,\cdots ,y_n)\approx \frac{1}{N}\sum _{j=1}^N{\hat{y}}_{n+i}({\theta }^{(j)}).$$
2. Use the fact that ${\displaystyle f_{\theta |y_1,\cdots ,y_n}(\theta )}$ is usually highly concentrated around ${\displaystyle \hat{\theta }=(\hat{\beta },\hat{\sigma })}$. We then ignore the small uncertainty of ${\displaystyle \theta }$ around ${\displaystyle \hat{\theta }}$: $$\mathbb{E}(y_{n+i}|y_1,\cdots ,y_n)\approx {\hat{y}}_{n+i}(\hat{\theta }).$$